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Abstract 

We study the dynamical response of a two-dimensional Ising model subject to a square-wave 
oscillating external field. In contrast to earlier studies, the system evolves under a so-called soft 
Glauber dynamic [P. A. Rikvold and M. Kolesik, J. Phys. A: Math. Gen. 35, L117 (2002)], for 
which both nucleation and interface propagation are slower and the interfaces smoother than for 
the standard Glauber dynamic. We choose the temperature and magnitude of the external field 
such that the metastable decay of the system following field reversal occurs through nucleation 
and growth of many droplets of the stable phase, i.e., the multidroplet regime. Using kinetic 
Monte Carlo simulations, we find that the system undergoes a nonequilibrium phase transition, in 
which the symmetry-broken dynamic phase corresponds to an asymmetric stationary limit cycle 
for the time-dependent magnetization. The critical point is located where the half-period of the 
external field is approximately equal to the metastable lifetime of the system. We employ finite-size 
scaling analysis to investigate the characteristics of this dynamical phase transition. The critical 
exponents and the fixed-point value of the fourth-order cumulant are found to be consistent with 
the universality class of the two-dimensional equilibrium Ising model. As this universality class 
has previously been established for the same nonequilibrium model evolving under the standard 
Glauber dynamic, our results indicate that this far-from-equilibrium phase transition is universal 
with respect to the choice of the stochastic dynamics. 
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I. INTRODUCTION 



Kinetic Ising or lattice-gas models with stochastic dynamics have been successfully ap- 
plied to study a number of dynamical physical phenomena, including metastable decay 
[l|, H, H, 3, H, JgIL hysteretic responses 0, H, 0], and magnetization switching in nanoscale 
ferromagnets [lO|, lul- Among the dynamic phenomena in such models that have attracted 
particular attention in recent years, is the dynamic phase transition (DPT) observed in 
systems with Ising-like symmetry that are driven far from equilibrium by an oscillatory 
applied force (typically a magnetic field or (electro) chemical potential). This phenomenon 

and later studied 
and 



was first observed in kinetic simulations of a mean-field model |12l. 
intensively by mean-field [11, H H EJ, Monte Carlo 0, S, UM, 



13| 



18, QJ, 20 



analytical |2^, |23|, |2J] methods. In this transition, the dynamic order parameter, which 



is the cycle- averaged magnetization, vanishes in a singular fashion at a critical value of the 
period of the applied field. Recently, strong experimental evidence has emerged that this 



nonequilibrium phase transition is observable in magnetic thin-film systems j25|, and an 
analogous phenomenon has been observed in simulations of a model of the heterogeneous 
catalytic oxidation of CO [2| 
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Perhaps the most fascinating aspect of this far-from-equilibrium phase transition is that 
it belongs to the same universality class as the corresponding equilibrium Ising model. This 
result is predicted from symmetry arguments jls], 2^ and has been confirmed by exhaustive 
kinetic Monte Carlo simulations 0, S, @, [HI and analytical results [21,]. Very recently, the 
field conjugate to the dynamic order parameter was identified as the cycle-averaged applied 
field, and a fluctuation-dissipation relation valid near the nonequilibrium critical point was 
numerically established 

The physics of equilibrium phase transitions is well understood, and it is well estab- 
lished that structures arising from different dynamics that obey detailed balance and respect 
the same conservation laws exhibit universal asymptotic large-scale features. However, the 
mechanisms behind nonequilibrium phase transitions are not that well known, and the de- 
pendence on the specific dynamic is still an open question. Except for the study of the 
model of CO oxidation [2^, 23], all previous kinetic Monte Carlo simulations in which this 
DPT was observed, were performed with the standard stochastic Glauber [30| or Metropolis 
'm dynamics. All these studies, including the study of CO oxidation which used a very dif- 



ferent dynamic, found critical exponent ratios consistent with the equilibrium Ising values, 
j/u = 7/4 and jS/u = 1/8. This gives a reasonable indication that this DPT is universal 
with respect to details of the model and the stochastic dynamics. However, a more direct 
test of just the universality with respect to the dynamics would be to use a significantly 
different stochastic dynamic for the two-dimensional kinetic Ising model. Such a test is the 
subject of the present paper. 

All stochastic dynamics that respect detailed balance eventually lead to thermodynamic 
equilibrium [32|, and all dynamics that obey the same conservation laws also give the same 
long-time dynamics (e.g., a t^^^ dependence of the characteristic length for phase ordering 
with a non-conserved order parameter and a t^^^ dependence for phase separation with a 
conserved order parameter [3^1]). However, it has recently been demonstrated that different 
stochastic dynamics give quantitatively dramatically different results for low-temperature 
nucleation 

la 



33, 1381, m 



35|. as well as for the nanostructure and mobility of field-driven interfaces 



40 



as 



The differences are particularly striking between dynamics known 
hard," in which the effects of the configurational and field-related ("Zeeman") energy 
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contributions in the transition rate do not factorize, and "soft," for which such factorization 
is possible (s^ . [sgI . 4i| . The class of hard dynamics includes the standard Glauber and 



Metropolis dynamics, while the soft dynamics here will be represented by the "soft Glauber 



dynamic" introduced in Ref. 36|], whose transition rate is given in Sec. [TTl Briefly, the 
nanostructure of field-driven "hard" interfaces is characterized by a local interface width 
and mobility that increase dramatically with the strength of the applied field, while "soft" 
interfaces remain relatively smooth and slow- moving, independent of the field [SQ]. Similarly, 
low-temperature nucleation under hard dynamics becomes very fast for strong fields, while 
under soft dynamics it remains thermally activated and thus very slow, even for very strong 
fields js^]. While the soft Glauber dynamic is probably not particularly relevant to any 
specific physical system, it is ideally suited for comparison with the standard, hard Glauber 
dynamic in investigating universal properties of the DPT. 

The rest of this paper is organized as follows. The kinetic Ising model and its dynamics are 
introduced in Sec. [TTl and the numerical results and finite-size scaling analysis are presented 
in in Sec. IIIII A summary and conclusions are given in Sec. IIVI 



II. MODEL AND DYNAMICS 

For this study we choose a kinetic, nearest-neighbor, Ising ferromagnet on a square lattice 
with periodic boundary conditions. The Hamiltonian is given by 

n=-JJ2^^^:i-H{t)J2'i^ (1) 

where Sj = ±1 is the state of the spin at the site z, J > is the ferromagnetic interaction, 
"^(^ijx runs over all nearest-neighbor pairs, runs over all lattice sites, and H(t) is an 
oscillating, spatially uniform applied field. The magnetization per site 

1=1 

is the density conjugate to H(t). The temperature T (in this paper given in units such that 
Boltzmann's constant equals unity) is fixed below its critical value {ksTc = J/ ln(l + "\/2)), 
so that, when there is no external field, the system has two degenerate equilibrium states 
with magnetizations of equal magnitude and opposite direction. When an external field is 
applied the degeneracy is lifted, and the equilibrium state is the one with magnetization in 
the same direction as the field. If the external field is not too strong, the state with opposite 
magnetization direction is metastable and eventually decays toward equilibrium [1]. This 
model is equivalent to a lattice-gas model with local occupation variables Cj = (sj + l)/2 G 
{0, 1} and (electro) chemical potential n oc H (for further details, see Ref. |37l]). 

The system evolves under the soft Glauber single-spin-fiip (non-conservative) stochastic 
dynamic with updates at randomly chosen sites. In the lattice-gas representation, this 
corresponds to adsorption/desorption without lateral diffusion. The time unit is one Monte 
Carlo step per spin (MCSS). When the system is in contact with a heat bath at a temperature 
T, each proposed spin flip is accepted with probabihty 

W^c = . (3) 
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Here /3 = 1/T, /S.Ej is the energy change corresponding to the interaction term, and /S.Eh 
is the energy change corresponding to the field term in the Hamiltonian, Eq. ([1]). This 
transition probabihty is to be contrasted with those of the standard, hard Glauber dynamic. 



= T-T7-^ , (4) 

1 + exp(/yA_£/) 

and the Metropolis djTiamic, 

Wm = Min[l, exp{-(3AE)] , (5) 

where AE = AEj + AEh is the total energy change that would result from a transition. 

The dynamical order parameter is the time-averaged magnetization over the kth cycle of 
the oscillating field [121 ], 



Qk = 7;— dtm{t) , (6) 

where ti/2 is the half-period of the applied field. The cycle is chosen such that it starts when 
H(t) changes sign. We also measured the normalized cycle-averaged internal energy, 

E 1 '■'=2* 



J(k-l)2Uio ^ /.-x 



J 2ti/2 y(fc-i)2ti/2 " 

As previous studies indicate 0, 0, 0] , the DPT transition essentially depends on the com- 
petition between two time scales: the average lifetime of the metastable phase, {t{T, Hq)), 
and the half-period of the applied field, ti/2- The metastable lifetime (r) is defined as the 
average time it takes the system to leave one of its two degenerate zero-field equilibrium 
states, when a field of magnitude Hq opposite to the initial magnetization is applied. In 
practice the metastable lifetime is measured as the first-passage time to zero magnetization. 

It is well known that metastable Ising models decay by different mechanisms depend- 
ing on the magnitude of the applied field Hq, the system size L, and the temperature T. 
Detailed discussions of these different decays regimes are found in Ref. More recently 
it has also been shown that, contrary to some common beliefs, there is also a strong de- 



pendence on the specific stochastic dynamics [3J, l35|. For the purpose of this study the 
temperature, the system sizes and Hq are chosen such that the metastable phase decays by 
random homogeneous nucleation of many critical droplets of the stable phase, which grow 
and coalesce, the so-called multidroplet (MD) regime. The metastable lifetime in the MD 
regime is independent of the system size 



III. MONTE CARLO SIMULATIONS 



The numerical simulations reported in this work are performed on square lattices with 
L between 64 and 256 at T = 0.8Tc. The system is subjected to an square-wave field H{t) 
of amplitude Hq = 0.3 J. The metastable lifetime was measured to be r = 145 ± 1 MCSS, 
almost twice as long as for the kinetic Ising model evol ving according to the standard (hard) 
Glauber dynamics under the same conditions 0, M, l20| . This is consistent with earlier 



observations of slow nucleation [34|] and interface growth [38j with this dynamic 
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The system was initialized with all the spins up, and the square-wave external field started 
in the half-period in which H = —Hq. After the system relaxed, the magnetization and 
energy reached a limit cycle (except for thermal fiuctuations), and all the period-averaged 
quantities became stationary stochastic processes. We discarded the first 2000 periods of 
the time series to exclude transients from the stationary-state averages. 

The time evolution of the magnetization is shown in Fig. [TJ For slowly varying fields 
(Fig. [T](a)), the magnetization follows the field, switching every half-period. In this region, 
Q ~ 0. For rapidly varying fields (Fig. [T](b)), the magnetization does not have time to switch 
during a single half-period and remains nearly constant for many successive field cycles. As 
a result, the probability distribution of Q becomes bimodal with two sharp peaks near the 
system's spontaneous equilibrium magnetization, ibmsp(T), corresponding to the broken 
symmetry of the hysteresis loop. The transition between these two regimes is characterized 
by large fiuctuations in Q. This behavior of the time series Qk-, shown in Fig. [21 is a clear 
indication of the existence of a dynamical phase transition between a disordered dynamic 
phase (the region where Q ~ 0), and an ordered dynamic phase (where Q 7^ 0). Notice that 
the transition occurs at a critical value Be = ^1/2/''" ^^^^ very close to unity, the value at 
which the half-period of the external field is equal to the metastable lifetime of the system. 
To further explore the nature of the DPT, we perform a finite-size scaling analysis of the 
simulation data. 



A. Finite-size scaling 

Previous studies indicate that although scaling laws and finite-size scaling are tools de- 
signed for equilibrium systems with a well known Hamiltonian, they can be successfully 
applied to far-from-equilibrium systems like the one we are analyzing here 0, i, i, H 



Since for finite systems in the dynamically ordered phase the probability distribution of 
the order parameter is bimodal, in order to capture symmetry breaking, the order parameter 
is better defined as the average norm of Q, i.e., To characterize and quantify this 

transition by using finite-size scaling we must define quantities analogous to the susceptibility 
with respect to the field conjugate to the order parameter in equihbrium systems. The scaled 
variance of the dynamic order parameter, 

xi = L\mL -m)i)- (8) 

has long been used as a proxy for the non-equilibrium susceptibility. A fiuctuation- 
dissipation relation was recently demonstrated, which justifies this practice by connecting 

to the susceptibility with respect to an applied bias field for a two-dimensional kinetic 
Ising model evolving under the standard Glauber dynamic [20| . 

In Fig. [3] we present the finite-size behavior of the order parameter and its fiuctuations. 
Fig. [3t^a) shows that this dynamic order parameter goes from unity to zero as ti/2 increases, 
showing a sharp change around t^^2 5 characterized by the peak in shown in Fig. [3](b). 
The absence of finite-size effects below the critical point is the signature of the existence of 
a divergent length scale. The height and the location of the maximum in xl change with L. 

In Fig. m we show the normalized time-autocorrelation function of the order parameter, 
defined as 

Q _ {Q(^)Q(^+p))-{Q(^))^ 
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The increasing correlation times with increasing system sizes are evidence of the critical 
slowing down of the system, and provides further support for the existence of a DPT. 
We also measured the period-averaged internal energy, Eq. ([7]), and its scaled variance 

xf = L'{{E'),-{E)l). (10) 

Both quantities are shown in Fig. [51 Again, in the absence of a fluctuation-dissipation 
relation, we use the scaled variance as a proxy for the analog of the equilibrium heat capacity. 
The correlation time was used to estimate the proper sampling interval for estimating the 
fluctuation measures and their error bars as described in Ref. [3^. 

It is very difficult to locate with precision the maxima of Xl Xl ^'^^ the individual 
finite system sizes. A more accurate estimation of the critical point at which the transition 
occurs in an infinite system can be obtained from the fourth-order cumulant intersection 
method. In Fig. [H]we plot the fourth-order cumulant Ul defined as 13 



as a function of ti/2 for several system sizes. Our estimate is t^^^ = (145 ± 1) MCSS, with a 
fixed-point value U* = 0.606 ± 0.004 for the cumulant (Fig. EJ^b)). The latter is consistent 
with the universal value for the two-dimensional equilibrium Ising model, U* = 0.6106901(5) 

Finite size-scaling theory for equilibrium systems (isl . 43 | predicts the following scaling 



forms at the critical point, 

(|Q|)^ocL-/^/^ (12) 

X? oc (13) 

which are also applicable to the far-from-equilibrium DPT 0, S, @, 20]. If the specific- heat 



critical exponent a = 0, as it is for the equilibrium Ising universality class, then we also 
expect the logarithmic divergence, 

xf^A + BHL). (14) 

These relations enable us to estimate the critical exponent ratios P/u and 7/z/ and verify 
the logarithmic divergence in the period-averaged internal energy fiuctuations. In Fig. [7] we 
present the results obtained by plotting the logarithm of {\Q\)l (Fig- El(a)), the logarithm 
of Xl (Fig- l^b)), and xf (Fig- E^c)); in term of the logarithm of L at ty2- We also plot 
the peak of the fiuctuations, x?(Peak) (Fig. EKb)), and xf (Peak) (Fig. [T^c)), since they 
asymptotically should follow the same scaling laws. After fitting the data with a weighted, 
linear least-squares algorithm, our estimates for the critical exponents are: (3/v = 1.44±0.06, 
7/z/ = 1.77 ± 0.04 (from the data at t^^2 )> ^^e data from x^(PeaA;) gives ^jv = 1.79 ± .02 
which agree within statistical error. Also, the straight line in Fig. [Tl^c) gives evidence of 
the logarithmic divergence of xf at the critical point. These results, together with our 
estimate for ?7*, give strong support to the hypothesis that the DPT observed is in the same 
universality class of the equilibrium two-dimensional Ising model. 
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IV. DISCUSSION AND CONCLUSIONS 



In this paper we have studied the dynamical response of a two-dimensional Ising model 
exposed to a square-wave oscillating external field. The system evolves under the so-called 
soft Glauber dynamic. In previous works it was established that, in the field and tempera- 
ture regions when the metastable decay occurs via a multidroplet mechanism, this system 
evolving under a standard (hard) Glauber dynamics undergoes a continuous phase transi- 
tion, with critical exponent ratios consistent with the equilibrium Ising values. The aim 
of the present study was to explore the universality of this far-from-equilibrium DPT with 
respect to the dynamics chosen to evolve the system. 

Our numerical results clearly indicate the existence of a DPT in the multidroplet regime. 
The transition depends on the competition between two time scales: the half-period of the 
applied field and the metastable hfetime of the system. We found that the metastable 
lifetime of the system evolving under the soft Glauber dynamics is roughly twice that of 
the same system evolving under the standard Glauber dynamic. However, in both cases the 
transition occurs at a critical point where both times are of the same order of magnitude. 
If the half-period of the applied field increases much beyond the metastable lifetime, the 
system is in a dynamically disordered phase characterized by a vanishing dynamic order 
parameter. A study of the autocorrelation function of the order parameter at the critical 
point provides evidence of critical slowing down, showing increasing correlation times with 
increasing system sizes. 

We applied the machinery of finite-size scaling, originally developed for equilibrium phase 
transitions, to estimate the critical point and the critical exponent ratios 13/v and 7/i/ 
for system sizes between 64 and 256 at T = 0.8Tc and Hq = 0.3J. Our estimates are 
a; 1.44 ± 0.06 and 7/1/ ~ 1.77 ± 0.04. These values are close to those of the two- 
dimensional equilibrium Ising model: 13/v = 1/8 = 0.125, 7/1/ = 7/4 = 1.75. Furthermore, 
our data strongly indicate a slow logarithmic divergence with L of the period-averaged energy 
fluctuations, consistent with the equilibrium Ising exponent a = 0. The flxed-point value of 
the fourth-order cumulant, U* , is also close to its expected universal Ising value, near 0.611. 

This study provides further evidence of the universality class of the dynamic phase tran- 
sition in kinetic Ising systems driven by an oscillating field, extending its domain to systems 
that evolve under different stochastic dynamics that lead to interfaces with significantly 
different structures on the nanoscale. 
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FIG. 1: (Color online) Time series of the magnetization (solid curves) in the presence of a square- 
wave external field (dashed lines), for two values of the half-period ti/2- ^1/2 = 400 MCSS, 
corresponding to a dynamically disordered phase, (b) ti/2 = 40 MCSS, corresponding to a dynam- 
ically ordered phase. The data were obtained for a system of size L = 128 at T = 0.8Tc and field 
amplitude Hq = 0.3J. 
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FIG. 2: (Color online) Time series of the order parameter Q for a system of size L = 128 at 
T = 0.8Tc and Hq = 0.3J. The horizontal trace near Q = +1 corresponds to a half-period of 
the field, ti/2 = 40 MCSS, well into the dynamically ordered phase (see Fig. mb)). The strongly 
fluctuating trace corresponds to ti/2 = 145 MCSS, very close to the DPT. The horizontal trace 
near Q = corresponds to ti/2 = 400 MCSS, well into the dynamically disordered phase. (See 
Fig. da)). 
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FIG. 3: (Color online) Dependence on the half-period ti/2 of the order parameter (a), and 

of its scaled variance xi (^)' shown for various system sizes, L. All the results correspond to 
T = O.STc and Hq = 0.3J. 
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FIG. 4: (Color online) Normalized autocorrelation function for the order parameter Q for ti^2 — 
145 MCSS at T = O.ST^ and Hq = 0.3J, shown for different values of L. 
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FIG. 5: (Color online) Dependence on the half-period ti/2 of the period-averaged internal energy 
(E) (a), and its scaled variance xf, (b) various system sizes, L. All the results correspond to 
T = O.STc and Hq = 0.3J. 
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FIG. 6: (Color online) (a) Dependence of the fourth-order cumulant Ul on the half-period ti/2, 
shown for various system sizes, L. (b) Enlargement of the region around the cumulant crossing. 
The horizontal and vertical dashed lines indicate the fixed point value U* ~ 0.606 and the critical 
half-period, t^^ = 145 MCSS, respectively. All the results correspond to T = O.STc and Hq = 0.3J. 
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FIG. 7: (Color online) Critical exponent estimates from the scaling relations. The symbols rep- 
resent the MC data, the straight lines are weighted least-square fits, (a) Calculating [3/u from 
Eq. (jl2p calculating 7/z/ from Eq. (jl3p . (c) The logarithmic divergence of the period-averaged 
energy fluctuations, based on Eq. ([E]). From data at t^^2 = 145 MCSS, circles, at the peaks, 
squares. 
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